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SUMMARY 

A simple  but  effective  system  for  the  manipulation  of 
polynomials  of  several  variables  in  APL  is  presented.  The 
system  is  especially  advantageous  in  situations  where  more 
sophisticated  symbolic  computing  systems  are  not  available, 
or  have  failed  to  solve  particular  problems.  The  system  is 
shown  to  successfully  solve  a problem  not  resolved  by  a more 
sophisticated  system. 
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1.  INTRODUCTION 

A quarter  century  of  development  in  the  area  of  symbolic  computing  has 
resulted  in  a wide  ranging  heritage  of  symbolic  computing  systems.  Excellent 
state  of  the  art"  reviews  of  this  area  have  resulted  from  the  1966  and  1971 
Symposia  on  Symbolic  and  Algebraic  Manipulation (ref . 1,2) . Considering  the 
unquestioned  power  of  these  systems,  it  is  perhaps  surprising  that  they  have 
not  gained  a wider  acceptance  within  the  scientific  community  as  a standard 
tool.  Part  of  the  problem  may  be  the  slow  process  of  education  as  to  what  is 
available.  The  categorisation  by  Moses (ref. 3)  of  systems  into  "conservative", 
liberal",  ''new  left",  and  "catholic"  eloquently  highlights  another  problem  - 
that  of  deciding  what  the  user  wants,  and  indeed  what  can  be  reasonably 
provided.  Systems  which  attempt  to  be  all  things  to  all  people  may  become  so 
complex  that  their  larger  computing  requirements  (time,  memory,  user  education 
and  systems  support)  discourage  their  incorporation  within  many  computing 
environments . 

There  is  therefore  a requirement  for  a smaller  scale,  more  "simple  minded" 
formalism  which  will  allow  isolated  researchers  to  harness  some  of  the  power  of 
symbolic  computing  even  though  they  may  not  have  access  to  for  may  not  know 

about)  the  more  extensive  systems  such  as  FORMAC,  ALTRAN,  REDUCE,  SCHOONSCHIP 
and  SIN.  ’ 

This  paper  presents  techniques  by  which  manipulation  of  polynomials  of  an 
arbitrary  number  of  variables  may  be  carried  out  using  APL.  A minimum  of  two 
simple  APL  routines  are  required  Cfor  polynomial  addition  and  multiplication) , 
but  the  simplicity  of  the  formalism  is  such  that  it  is  easy  to  extend  the  range 
of  functions  to  allow  for,  say,  polynomial  division,  and  the  manipulation  of 
rational  polynomials. 

APL  has  not  been  fully  exploited  as  a medium  for  symbolic  computing,  although 
its  potential  is  undoubted.  The  use  of  APL  for  manipulation  of  polynomials  of 
a single  variable  has  long  been  recognised,  and  it  is  now  included  in  basic 
texts(ref.4) . Other  users (ref. 5, 6, 7)  have  pointed  to  the  possibilities,  while 
generally  limiting  their  attention  to  specific  areas.  More  detailed  comments 
will  be  given  later  on  Surkan's  application (ref .7) . 

APL  is  especially  suited  to  symbolic  computing  for  several  reasons.  Firstly, 
the  programmer  does  not  have  to  worry  about  storage  control.  Routines  may  be 
coded  with  complete  generality,  without  regard  to  the  rank  or  size  of  arrays. 

The  storage  problems  caused  by  "intermediate  expression  swell"  in  some  other" 
systems  (see  Tobey(ref . 8))  and  'garbage  collection"  are  therefore  not  apparent  to 
the  user.  Secondly,  APL  functions  can  be  used  as  operators  which,  for  our 
purposes,  allow  polynomial  expressions  to  be  built  up  naturally  and  executed 
without  the  use  of  an  intermediate  interpreter.  Further,  APL  is  designed  to  be 
used  interactively.  This  is  a decided  advantage  when  the  form  of  a problem 
solution  is  not  known  in  advance.  The  user  is  in  a position  to  massage  the 
output  of  the  symbolic  computation  in  a way  which  he  determines  at  that  time. 

The  approach  and  the  routines  presented  in  this  paper  were  developed  to  solve 
a particular  problem  in  filtering  theory.  A solution  had  been  attempted  using 
FORMAC,  but  failed  because  an  unknown  intermediate  term  went  outside  the 
allowable  range.  This  raises  a further  advantage  of  this  simpler  approach. 

The  system  presented  here  is  sufficiently  transparent  so  that  if  something 
does  go  wrong  in  the  computations,  the  user  is  in  a good  position  to  work  out 
how  he  can  circumvent  the  problem  - either  by  reformulating  it,  or  by  writing 
his  own  routines  to  handle  it.  The  presentation  here  is  therefore 
deliberately  open-ended,  to  point  the  reader  towards  the  possibilities  rather 
than  to  spell  out  all  the  details. 

The  paper  assumes  a knowledge  of  APL. 
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2. 

Consider  a polynomial  in 


POLYNOMIAL  REPRESENTATION 
the  r variables  x^,...,xr. 

may  be  represented  as  an 

A = > • * • , ir) 

xi 


The  polynomial 
array  of  rank  r 

CD 


Thus  the  index  of  each  term  in  A is  1 greater  than  the  power  to  which  each 

variable  x.  is  raised.  Alternatively,  if  the  APL  index  origin  is  set  to  zero 
i 

(□  10-*-  0) , the  term  A[  i^ ; . . . ; i^]  gives  the  coefficient  of  the  term  in 
i i 

1 ...  x r.  This  representation  is  in  general  more  expansive  than  is 

strictly  necessary,  since  a large  proportion  of  the  elements  of  A may  be  zero. 

A more  compact  form  would  be  one  such  as  that  adopted  by  ALTRAN,  which,  for  N 
non  zero  terms,  would  represent  A by  only  N.(r+D  items (ref. 9) . This  more 
compact  form  was  also  adopted  by  SurkanCref .7)  in  his  APL  implementation. 
However  our  simplicity  of  representation  means  that  coding  of  routines  to 
manipulate  the  arrays  is  much  simpler. 

The  polynomial  form  CD  has  been  used  for  the  basic  arithmetic  operations 
add,  multiply  and  divide  between  polynomials,  but  it  is  capable  of  extension  to 
allow  manipulation  of  rational  polynomials,  using  as  building  blocks  the 
routines  to  handle  the  form  (1) . 

One  possible  representation  is  to  place  the  numerator  and  denominator  "back 
to  back".  Thus  the  numerator  is  represented  by  the  array 


A[  0, ,....]  , 


C2) 


and  the  denominator  by 


Mi;; ] 


Alternatively,  the  denominator  may  be  kept  as  a series  of  factors,  so  that 
the  numerator  is  represented  by 


A[  0; ; ] 


(3) 


as  in  C2) , but 


Mi;;...l,  i > 0 


are  factors  of  the  denominator. 

The  form  C3)  would  be  less  efficient  with  storage,  since  the  arrays 
Mi;;...]  would  all  have  to  be  filled  out  to  the  same  shape.  However,  in 
problems  where  the  same  factors  occur  frequently  in  denominators,  this  may  be 
less  demanding  than  searching  to  eliminate  common  factors  from  the  form  C2) • 


- 3 - 
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3.  POLYNOMIAL  MULTIPLICATION 

If  PA  and  PB  are  polynomials  represented  by  the  arrays  of  coefficients  A and 

B respectively,  then  the  polynomial  product  P^P^  can  be  found  by  forming  the 

outer  product  A°.xB,  and  summing  those  elements  of  the  outer  product  which 
correspond  to  the  same  exponents  in  the  polynomial  product.  The  procedure  for 
selecting  and  summing  these  elements  can  be  most  easily  visualised  as  a 
generalization  of  the  single  variable  polynomial  product. 

If  A and  B represent  polynomials  of  degree  n and  m respectively,  in  x,  the 
outer  product  will  be  an  (n+l)x(m+l)  array  of  products  (a.b.),  i = 0,...,n, 

j = 0,...m,  (where  A = (a  ....,aj  and  B = (b  , . . . ,b 
^ o n o m 


The  coefficient  of 


x in  the  polynomial  product  will  be  the  sum  of  all  terms  a.b.  such  that 

1 1 

i+j=k.  Thus  we  must  sum  all  elements  in  the  outer  product  along  the  reverse 
diagonal  i+j=k  for  each  k = 0,1,..., n+m  (see  figure  1). 

There  are  two  ways  in  which  the  elements  along  each  reverse  diagonal  of  the 
outer  product  may  be  selected.  The  first  is  to  generate  a selector  matrix 
which  has  the  same  shape  as  the  outer  product,  and  which  is  zero  everywhere 
except  for  l's  along  the  appropriate  reverse  diagonal.  The  APL  product  of  this 
selector  matrix  with  the  outer  product  will,  when  reduced  by  summation,  give  the 
desired  coefficient  of  a term  in  the  polynomial  product.  This  technique  has 
been  implemented  in  APL  and  has  been  found  to  be  adequate,  but  slower  than  would 
be  desirable.  The  slowness  is  due  to  the  need  to  separately  shape  the  selector 
matrix  for  each  term  in  the  polynomial  product. 

The  second  technique  which  has  been  implemented  is  performed  by  rotating  the 
outer  product  matrix  so  that  the  diagonals  i+j=k  become  columns,  which  may  then 
be  reduced  by  summation  in  a single  APL  operation. 

The  process  is  illustrated  for  the  single  variable  case  in  figure  2,  where 
the  outer  product  has  been  padded  with  n columns  of  zeros  to  the  left,  and  the 
rows  then  rotated  by  amounts  ranging  from  n for  the  first  row,  to  0 for  the  last 
row.  Reduction  by  summation  over  the  first  axis  then  gives  a vector  which  is 
the  array  representation  of  the  polynomial  product. 

For  the  generalization  to  polynomials  of  r variables,  assume  P and  P are  of 
a A B 

degree  n^,...^  and  in  each  of  the  variables  respectively.  The 

outer  product  A°.xB  will  have  shape  (n1+l)  , . . . , (n^l)  , (rn^l)  , . . . , (mr+l)  , and  the 

product  will  be  of  degree  = n.+nn,^! , . . . ,r.  If  we  pad  the  outer  product  in 

the  same  way  as  for  the  single  variable  case,  the  resulting  array  will  have  shape 


> • • • » (nr+l) > (P^+l) > • • • j (Pr+1) » 


that  is,  we  have  padded  n layers  of  zeros  onto  the  left  side  of  the  (r+i)th 
axis  for  each  i. 

A separate  rotation  will  now  need  to  be  carried  out  for  each  variable.  For 
the  first  variable,  this  means  defining  an  array  C of  shape 


such  that 


OVl)  >•••»  Or+l)  > (P2+l)  >•••>  (pr+l) 

C[  0; . . . ;]  is  n, 

C| 1; — j]  is  n -1, 


C[  n 


;]  is  0. 


and 
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After  the  rotation  has  been  performed,  each  section  of  the  array 
perpendicular  to  the  (r+l)st  axis  will  be  characterised  by  (ii  +>ji ) = constant, 
where  ii  is  the  index  of  the  first  and  ji  is  the  index  of  the  (r+l)st  axis 
respectively.  Each  such  section  would  have  the  same  shape  as  C. 

If  this  rotation  operation  is  now  carried  out  with  respect  to  the  second 
variable  in  the  same  way,  the  section  of  shape 

(n1+l) , . . . , (nr+l) , Cp3+1) > • • • » (pr+1) 

which  is  perpendicular  to  the  (r+l)st  and  (r+2)nd  axes  will  be  characterised  by 
ii+ji  = constant  and  h+h  •=  constant. 

After  similar  rotations  for  each  of  the  r variables,  each  section 
perpendicular  to  the  (r+1) st, . . . , (2r)th  axes  will  have  shape  (n^+1)  , . . . , (n^+1) , 

and  will  be  characterised  by 


i +i  = constant  = k , s 
s s s 


1>. 


,T 


Thus  elements  contributing  to  each  term  in  the  product  are  located  in 
columns  defined  by  the  first  r axes  of  the  array,  and  the  array  representing  the 
polynomial  product  may  be  formed  by  reduction  by  summation  over  each  of  these 
first  r axes. 

While  this  technique  is  a little  more  difficult  to  visualize  than  that 
involving  the  generation  of  a selector  matrix,  it  is  shorter  (30  times  shorter 
for  multiplication  of  polynomials  of  degree  9 in  two  variables)  because  the  bulk 
of  the  operations  are  carried  out  once  for  each  axis  (variable) , rather  than  once 
for  each  term  in  the  product. 


4.  POLYNOMIAL  DIVISION 

The  system  used  for  polynomial  division,  while  it  thoroughly  exploits  APL 
techniques  and  therefore  may  look  quite  different,  is  identical  in  concept  to 
that  outlined  by  Collins (ref . 10)  for  the  PM  system. 

Consider  two  polynomials  A and  B in  r variables.  We  assume  A is  of  degree 

n in  x.  and  B is  of  degree  m.  in  x..  If  n.^m.  for  any  i we  know  B will  not 
ii  iiii 

divide  A and  hence  without  loss  of  generality  we  may  assume  nu,i=l, . . . ,r. 

We  start  by  defining  the  most  significant  term.  This  is  done  by  defining  a 
primary  ordering  of  the  non -zero  terms  according  to  the  exponent  of  Xi  , a 
secondary  ordering  according  to  the  exponent  of  X2  and  so  on.  In  APL  this  is 
implemented  by  recursively  removing  outer  layers  of  zeros,  until  the  array  is  at 
its  smallest  significant  size,  and  then  dropping  all  except  the  last  layer  of 
the  remnant  array  in  the  first  dimension  on  the  first  cycle,  the  second 
dimension  on  the  second  cycle,  and  so  on  to  the  r-th  dimension  on  the  r-th  cycle. 

The  most  significant  term  in  B can  then  be  divided  into  the  most  significant 
term  in  A;  the  result,  q say,  added  to  the  quotient  Q,  and  a new  A'=A-Bq  formed 
in  the  usual  manner  of  polynomial  long  division.  The  process  will  halt  when 
the  leading  term  in  B no  longer  divides  the  leading  term  in  A'.  A'  will  then  be 
an  array  which  has  only  zero  elements  for  (i^>  m^,^^  > • • • > ^r'>  * Note  that 

this  does  not  necessarily  mean  that  the  degree  of  Q in,  say  xs,  is  less  than  ms. 

In  fact  Q may  still  be  of  degree  n^  in  each  of  the  variables  x^,...,x^. 
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(E.g.  consider  A = x2 +y2 +xy  and  B = xy)  . In  2 dimensions  this  may  be 
illustrated  by  representing  the  remainder  A’  as 

A'  =(*«  ^ \ 

VAi  0 J 

where  Ao  has  shape  n,m  and  Ai  and  A2  are  not  necessarily  zero. 

5.  FURTHER  CONSIDERATIONS 

The  basic  arithmetic  operations,  multiplication,  division  and  addition  are 
sufficient  to  solve  a wide  range  of  polynomial  manipulation  problems.  Listings 
of  APL  routines  to  perform  these  functions  are  given  in  Appendix  I. 

The  most  expensive  of  the  basic  routines  ADD,  BY,  and  DIV  in  terms  of  storage 
requirements  is  the  polynomial  multiplication  routine,  BY.  This  routine  must 
accommodate  an  array  equal  in  size  to  the  outer  product  of  one  factor  and  the 
polynomial  product.  Under  normal  circumstances,  and  especially  in  a virtual 
storage  environment,  this  should  not  cause  problems.  However  the  applicability 
of  the  routine  can  be  extended  if  a checkpoint  is  added  before  the  formation  of 
the  outer  product  in  BY,  at  which  the  available  working  area,  □ WA,  is  tested 
against  requirements.  If  necessary  one  of  the  factors  can  then  be  split 
(presumably  half  way  along  its  longest  axis)  for  a recursive  call  to  BY.  If 

this  occurs  then  the  section  of  BY  which  processes  outer  products  will  be 
bypassed  in  the  primary  call.  Because  lowest  powered  exponents  are  eliminated 
in  BY  before  formation  of  the  outer  product,  this  procedure  may  considerably 
reduce  the  workload  if  either  of  the  factors  is  large  and  has  widely  spaced 
non-zero  elements.  The  listing  of  BY  in  Appendix  I includes  such  a checkpoint 
with  a rudimentary  test  of  available  space.  Tests  carried  out  on  the  self- 
multiplication of  a 10  x 10  non-zero  array  showed  that  the  execution  time  was 
increased  by  about  50%  when  recursion  was  carried  to  the  fourth  level. 

There  are  further  operations  which  may  usefully  be  performed  by  using  the 
basic  routines  as  building  blocks. 

One  desirable  facility  is  to  be  able  to  eliminate  common  factors  from  two 
polynomials,  say,  numerator  and  denominator.  A traditional  approach  to  this 
has  been  the  search  for  a greatest  common  divisor.  In  theory  this  can  be 
achieved  by  a straightforward  application  of  Euclid's  theorem,  but  in  practice 
it  has  proved  to  be  one  of  the  great  problem  areas  of  symbolic  computing  (see, 
for  example,  Collins(ref. 11)) . 

A simpler  but  more  limited  approach  that  has  been  used  successfully  by  the 
author  has  been  to  build  a "library"  of  factors  which  arise  naturally  in 
denominators  in  the  course  of  a series  of  calculations,  and  to  look  for  factors 
from  amongst  these.  If  a series  of  routines  (say  ADDR,  BYR  and  OVER)  are  built 
up  from  the  basic  ADD  and  BY  to  handle  rational  polynomial  functions,  then  it 
would  be  appropriate  in  OVER  to  check  to  see  if  a new  factor  could  be  added  to 
the  library  of  denominator  factors. 

Another  useful  facility  is  the  ability  to  substitute  new  expressions  for 
variables.  A routine,  REDRNK,  has  been  written  to  perform  this  and  its 
listing  is  included  in  the  appendix.  The  only  limitation  is  that  the  user  may 
need  to  transpose  co-ordinates  of  the  input  polynomial  so  that  substitution  is 
made  for  the  correct  variable  (the  last  one).  If  the  substituted  polynomial 
includes  a new  variable  the  input  polynomial  will  also  need  to  be  recast  so  that 
it  is  a function  (albeit  trivial)  of  that  variable.  Both  these  operations  are 
straightforward  in  APL. 

Further  functions  such  as  polynomial  differentiation  can  also  be  easily 
coded  in  APL. 
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6.  AN  APPLICATION 

The  following  application  is  one  which  arose  in  the  course  of  research  in  the 
area  of  Kalman  filtering  and  smoothing.  The  aim  was  to  determine  the  band- 
widths  of  a particular  Kalman  filter-smoother  as  a function  of  the  bandwidths  of 
the  Kalman  filter  without  the  smoother. 

Under  steady  state  conditions  the  filtering  and  smoothing  matrices  may  be 
described  in  terms  of  3 parameters  f , f and  f , which  are  proportional  to  the 

S V 3. 

filter's  initial  response  to  a step  in  position,  velocity  and  acceleration 
respectively.  In  fact,  f will  approximate  the  (non-dimensional)  frequency 

response  of  the  filter  with  respect  to  position.  (The  same  cannot  be  said  of 
fv  and  f because  of  phase  shifts.) 

Only  one  of  these  parameters  may  be  specified  independant ly,  and  we  may  write, 
for  t = V(TTS) , 

fs  = Ot2)  ; fv  = 2 (1-t)2  ; fa  = (1-t)3  /(1+t) 

If  bs,  b , b&  are  the  respective  smoothing  responses,  it  is  possible  to  show 

that  they  may  be  defined  so  that  they  are  the  diagonal  elements  of  a 3 x 3 matrix 

W,  given  by 

W = (PS-$)(I-A_1) 

where  $ = / 1 1 % \ A"1  = / a/d  -1  h \ 

Oil  b/d  1 -1 

\0  0 1 / \ 2c/ d 0 1/ 

f , 2f  \ 

v a \ 

f (12-10f  -f  )/(8f  (1-f  ))  , f (2f  -f  )/(l-f  ) 

v s vJ  y s v s"  a^  s vJ  K sJ 

f ((2f  -f  )/(l-f  ) , 2f  f /(1-f  ) / 

a'.v  s vj/  k sj  > a v ' K SJ  / 

a = (l-fv+fa);  b = (fv-2fa);  c = f&;  and  d = (l-fg) , 

and  S is  the  solution  of  the  equation 

S-ATSA  = r $ 

where  F = / 1 0 0 \ 

0 0 0 
\ 0 0 0/ 

Here  P is  proportional  to  the  covariance  of  the  filtered  state  vector,  F is  a 
scaled  information  matrix,  and  $ the  transition  matrix  for  a unit  data  rate. 

The  only  significant  remaining  problem  is  the  solution  of  the  equation  for  S. 
Since  A is  algebraically  simpler  than  A,  this  may  be  equivalently  expressed  as 

(A-1)1  S(A_1)  - S = (A-1  )T  T $ (A-1  ) 
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It  was  the  solution  of  this  set  of  9 equations  in  the  9 unknowns 
Sij ,i, j=l,2,3  which  was  attempted  using  the  standard  procedures  of  FORMAC  and 

which  failed  because  some  intermediate  term  had  gone  beyond  the  allowable 
magnitude  range.  Solution  of  the  equations  was  then  attempted  by  hand. 
Although  a range  of  relatively  simple  subrelations  could  be  derived,  they  were 
still  too  cumbersome  to  solve  by  hand. 

The  relations  were 


Si  1 = Si  2 ■ + S2 1 

Si  2'  = 2 (Si  3 + S2  2 ) 

S2 1 ■ = 2 (S3 1 + S2  2) 

8CS3  3 = 1 + 2(d-l)  (S13+S31)' 

S31  • = S3  2 = - S23 

2(l+d)  (S13.+S31  ) + 2 (2a+b+2d)  Sir  = 1 

2 (b+c-d(b+3c))  (Si  3 - S31)  = b+c 

and  2(b+c-d(b+3c))  S3j  ■ = -d 

Solution  was  again  attempted  using  FORMAC,  which  again  failed. 

The  subrelations  were  then  coded  in  APL  using  the  polynomial  manipulation 
routines  ADD  and  BY  and  the  results  were  used  to  calculate  the  responses  b 

bv>  ba5  as  Pairs  of  polynomials  (numerator  and  denominator)  in  two  variables  t, 

and  6 = (1+t)-1 . The  solutions  that  were  obtained  are  listed  in  Appendix  II. 

The  validity  of  the  polynomial  solutions  was  confirmed  by  the  numerical 
solution  of  the  equations  for  S and  W for  a range  of  fg.  However,  the  primary 

form  of  the  polynomial  solution  as  numerators  WN11,  WN22,  WN33,  and 
denominators  WD11,  WD22,  WD33,  was  still  too  complex  to  give  insight,  and  it  was 
therefore  desirable  to  simplify  these  expressions  by  use  of  polynomial 
manipulation. 

The  first  step  was  to  eliminate  common  factors  by  testing  each  numerator  and 
denominator  against  a library  of  9 factors  which  arose  in  the  course  of  the 
calculations.  The  second  step  was  to  make  use  of  the  rational  function  nature 
of  the  result  to  eliminate  the  variable  9.  This  was  done  first  by 
substituting  ? = 1/0  (effected  by  reversing  the  order  of  the  last  co-ordinate  in 
the  numerator  and  denominator)  and  then  by  substituting  f = (1+t) . After 
testing  the  results  for  factors  (1+t)  and  (1-t) , it  was  found  that 

bs  = (l+6t-34t3 +39t4 -4t5  -8t6)/(t2  (l+t2  ) (l+4t+t2  ) ) 

by  = (1-t2 )/ (l+4t+t2 ) 

and  ba  = 4t(l-t)/((l+t)(l+4t+t2)) 

and  it  is  remembered  that  t2  = 1-f  . 

s 

The  output  of  the  two  steps  in  the  simplification  is  included  in  Appendix  II. 


7.  DISCUSSION 

Some  points  of  comparison  should  be  made  between  the  APL  system  of 
polynomial  manipulation  presented  in  this  paper,  and  those  proposed  previously 
by  Surkan(ref . 7) . In  his  paper,  Surkan  reported  the  development  of  APL 
routines  to  add,  subtract,  multiply  and  differentiate  polynomials.  He  used  a 
different  form  of  polynomial  representation  in  which  each  term  is  represented  by 
its  coefficient  and  the  exponents  of  the  (ordered  set  of)  variables. 


ERL-0004-TR 


8 


Tests  were  run  with  Surkan's  example  of  orbit  calculation  to  compare  the 
effectiveness  of  the  two  approaches. 


N 


N 


es  U,  V, 

, and  W,  given 

that  Fo 

+ ^FN-l  + 

^N-l 

3U 

3V 

3W 

"Sj-i 

+ v55n-i  ♦ 

31) 

0V 

3w 

N-l 


N-l 


and  U = -3UV,  V = W-2V,  W = -V(U-2W) . (Note  that  the  equation  for  V was 
listed  incorrectly  in  Surkan's  paper,  although  it  was  programmed  correctly  in 
his  figure  2.) 

The  current  set  of  routines  performed  these  calculations  to  N = 19  in  18.2  s 
on  an  IBM  370/168  with  VSAPL.  This  compares  most  favourably  with  Surkan's 
quoted  time  of  20  min  on  an  APL  online  terminal  to  an  IBM  360/50,  in  spite  of 
the  difference  in  systems. 

t>  of  more  consequence  is  the  difference  in  the  results  at  N = 11.  Although 
some  terms  agree  with  those  of  Surkan's  figure  3,  most  do  not.  Only  the 
addition  routine  was  listed  in  Surkan's  figure  1,  and  so  it  is  difficult  to  be 
more  precise  about  where  the  difference  arises.  Since  the  same  equations 
programmed  in  FORMAC  give  answers  identical  with  those  from  the  current  set  of 
routines,  it  must  be  assumed  either  that  this  author  has  misinterpreted  the 
example,  or  else  the  routines  reported  in  reference  7 were  in  some  way  faulty. 

A listing  of  the  correct  solution  for  N = 11  is  given  in  Appendix  III. 

It  is  interesting  to  note  that  the  FORMAC  solution  took  5.65  s in  primary 
execution  time  and  8.43  s including  the  preprocessing,  compilation,  and  link- 
editing steps.  This  does  not  diminish  the  value  of  the  APL  approach  for  the 
reasons  discussed  in  section  1.  In  fact  the  ratio  of  execution  times,  at  a 
little  more  than  2:1,  appears  to  be  substantially  less  than  the  10  - 100:1 
generally  accepted  for  execution  times  of  problems  coded  in  APL  to  those  for 
coding  in  other  higher  level  languages  (such  as  FORTRAN  and  PL/1) (12) . The 
advantage  with  APL  normally  comes  with  the  time  to  program  and  debug,  and  this 
will  dominate  in  programs  which  do  not  require  repeated  execution.  It  is  t is 
'problem  solving'  orientation  of  APL  which  is  especially  suited  to  symbolic 

computing. 
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APPENDIX  I 

APL  PROGRAM  LISTINGS 


v/mmv 

7 RES*- A BY  P ; HIO  ;C  ;RM  ;I  ;J  }V  ; RNK  }RA  ;RB;FAC 

ri]  ,77^(Dp4)r(pp7)rnio^-o 

r 2 1 pi  CHECK  FOR  HULL  OP.  TRIVIAL  POLYNOMIAL 

f 31  ->((0  = + /,  U)v(o  = + /,  |B)  )/NULL 

[4]  -»((0<l /pP*-SmRIP  B)a(  0<l/pA*-S?RIn  A)  )/*?. 

[51  7 u L L : ->0  x p p p ES*-(  RNK  p 1 ) p 0 

[fi  1 TA  :-*-OxppREE*-Rx  + / ,A 
[7]  ?Ri-*>QxppPES*-A*  + / ,3 

m ^ i*r /p>n ,(i=r/p^) ,((op/i)=pPB) )/(?/, r/?,. 9i) 

[9]  B*-((  (RNK-ppB)pl)  ,pB)pB 

[101  A*-((  (RHK-ppA)pl)  ,pA)pA 

[11]  a REMOVE  COMMON  TONERS  FROM  A AND  3 

[12]  Si  :RB*-(pB)  + 0xRA*-(pA  ) +FAC*-R'J.Kp  I*-0 
[13  ] £0:  7-e-I+l  + i ( RNK-I+ 1 )+J«-0 

[14]  L4  : ->-(0  = + / , | (/M[i_n,(j4v7+1)  , ?M  t .V 1 ) +4  ) / r,4 
[15  1 «74-0x^[J>ff4m -FAC[Il*-J- 1 

[16]  4<-(-/M)+4 

[17]  £J?j+(0  = + /t  I (TPClIl  ,(JW  + 1)  ,/?P[B7])+5)/L5 
[IP]  PBl.Il+RBtTl-J+J-l 

rig]  P*-(-RB)+3 

[70]  -*■  ( RNK > I*- 1 + 1 + 0 x FA  C [ I ] *-FA  C [ J ] + J ) / L 0 

[21]  A CHECK  AVAILABLE  WORK  A^EA 

[22]  ->•(  0<H[M-40  + (2  0xmY)  + x/2  0,P/  ,RM*-RA+RB-1  )/52 
[231  I*t-l  + 0x«N-[/7M 

[24]  L:-y(J*(I-l)UI*-I+l)  + RA)/I 

[25  ] RES*-(  -U  + pP  ES  ) + RES  *-B  B v (N*-(  ( 7-1  ) p 0 ) ,(«7«-[  0 . 5x,7)  ,(RVK~I)o 
0 ) \A 

[26]  ->(  RNK-pRM*-o  RES*-PES  ADD  B BY  ((  (I -1)  t»M)  ,J  ,T  *PA)  +/  ) / RES" 

[2  7]  S2  :I*-0*ppRES<-(  *A  , ( -RM ) ) ♦/.  <>  ,xB 

[2  8]  L 1 : N*-  ((IAEA)  , 1 , ( ( 1+1  ) 4 ) , (7+W)  , (I+1)*PM) 

[29]  ( 7?>l  [T  ] = 1 Al\pC*-N  pj  *-PA  [ T]  -1  ) /EQ1 

[30]  12  : + (RAin>l+npC*-C  JIlNpI+J-l  )/L2 

[31]  EQ1  :-*■((  I *-I  + l)<0.5xpp  RES  *-C<k  [ 7?  NK +T 1 RES  ) /Li 

[32]  L 3 : -*-(  PNK < p p PES*-+  / [ 0 IRES  )/L  3 

[33]  a RESTORE  COMMON  POWERS  ^O  RESULT 

[34]  REST  : RES*-STRIP  ( -PM+FA  C ) t RES 

7 

7 /.  0 jO  [ f]  ] 7 
V RES*- A ADD  B \R 

[1]  a ENSURE  A AfW  B HAVE  THE  SAME  RANK 

[2]  ->(  (pp4)  = ppJ?)/71 

[3]  R*-  ( p p4  ) [ p p5 

[4]  ]-^((  (7-pp])ol)  ,p3)pl 

[5]  A*-(  ( (R-opA  )pl  ) ,pA  )pA 

[fi]  a PAD  0Um  SMALLER  OF  EACH  DIMENSION 
[ 7 ] 51  : RES*-STPIR(RAA  ) + ( *>*-(.  pA  ) [ pB  ) +7? 

7 

VSTRIPtWiR 

7 R-t-STRI13  A-,PA  ;NNK-,I-,J-,AUA2 

[1]  a REDUCES  A m0  SMALLEST  SI  ON IF.  SIZE 

[2]  -*-((+  / , U)=  + /,.fl«-((l[PJ7K*-pP-4  )pl)p4)/0 

[3]  T-«-l+e7-<-0xpp  R*rA 

[4]  Ll:Al*-(I-l)  + RA*-pR 

[5]  A2*-I  + RA 

[6  1 L 2 :-*■(<>  = + /,  I (4.1  ,( -JW  + 1 ) ,A2)+P)/L2 

[7]  -*-(  (Tf-J+1  ) <pp. *?*-((  (T-l  )pJ*-0  ) , (1-J) , ( RNK-I)pO)\D) /T,i 

7 
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VDIVl n]V 

V RES+A  DIV  B -,I  lUTO  ;?NKiQ  -,AS  iRS  ',CB  ;F?C3  iRDQiRSQ 

[1]  B ?0L. A i POL.B , RETURNS  QUOTIENT  AND  REMAINDER 

[2]  r!J3«-l+0x+/  ,RNK*-(p  pA+STviP  A ) fppB*-SmRIR  B 

[3]  -»-(  ( p pA  ) =p  pB  ) /SO 

[4]  4<-((  (RNK-ppA)pl)  ,pA)pA 

[5]  3«-((  (RNK-ppB)pl)  ,pB)pB 

[6]  50:-*-(l<r/pB)  /PI 

[7  ] ->0xpp/?ff5^-(2  ,p<?)  + (l  ,p<2)p<7«-(/U  + /,B  ) 

[8]  ft  FIND  THE  HIGHEST  POWERED  TERM  I”  3 

[9]  Sl:RCB+-RNKpO*I+-l 

[10]  RS+pCB+-B 

[11]  LB  :RCBin*-RCBlIl+RSin 

[12]  -*(  ( J-+-J+1  ) <pRS-+-pCB+-STRIP  CB+(((I-1)*RS)  ,~1  tI*RS)*CB) /LB 

[13]  0+(RNK pl)p0 

[14]  B CAN  A QUOTIENT  BE  FOUND? 

[15]  LI  :-K0>L/((  pA  )-pB)  )/NO 

[16]  +(  0 = + / , \AS+(RSQ+(RCB)  -I+-1)\A) /NO 

[17]  RS+pAS+STRIP  AS 

[18]  L 2 : RSO  [ I ] +-RSQ  [ I ] +PP  [ I ] 

[19]  ->(  (I+I  + l^pPP^-p/tS^tfjp  /?  P«-(  ( (J-l  ) + Pp  ) , ~1  ,J+p,?)+45)/P2 

[2  0]  b RSO  IS  HIGHEST  POWER , 4/7#  /IP  ITS  COEFFICIRN m 

[21]  RDQ+RSQ-RCB-l+I+O 
[2  2]  4P-*-+/ ,/lP*PB 

[23]  /IPP  P0<?p(  ((*/PP#)-l)pO),i4P 

[24]  +(.RNK=ppA+STRIP  A ADD-( 1 -RDO+pB ) \B*AS ) /LI 

[25]  NO:RES+(RS*Q)  , [ 0 . 5 ] ( PP<-(  p 0 ) T p/l  ) +/1 

V 

NRATPOLimv 

V PPP+/1  P/l  B ; □ JO  ; J ; P.V£  ; P/1  ;RB  ;JA  ; JR  ;N 

[1]  RNK+(ppA)r(ppB)n 

[2]  pi  ENSURE  A AND  B HAVE  THE  SAME  RANK 

[3]  -*(  (RNK=ppA  )*RNK=ppR  )/Sl 

[4]  B*(( (RNK-ppB)pl) ,pB)pB 

[5]  A*-(  ( (RNK-p  pA  )p  1 ) ,pA  )pA 

[6]  Pl:-*((0  = + /,  M)v(o  = + /,  \B))/S2 

[7]  b REMOVE  COMMON  POWERS  OF  VARIABLES 

[8]  I+DIO+-0X  + / (RA+pA)  ,RB+-pB 

[9]  Li  :N*-I+l  + \ (RNK-I+1 ) +JA+-JB+-0 

[10]  LA  :-*(0  = + / , | (RA  [i  J]  , (JA+-JA  + 1)  ,RA  [P])-M  ) /LA 

[11]  LB:+( 0=+/, | ( RB [ 1 1 ] , (JB+JB+1 ) ,RBZN ] ) +P ) /LB 

[12]  RA  [ J ] +-RA  [J]-e7/l-»-(t7/lLt7P)-l 

[13]  RB  [ J ] *-RB  [ J ] - J/l 

[14]  /l*-(-P/i)+/l 

[15]  B+(-RB)iB 

[16]  +(RNK*I+I+1 )/Ll 

[17]  P2  : RES-*-(  2 ,RA)p  ( ,RA  ) , , (P/l<-(  pA  ) TpP  ) +B 

V 

vpppppprmv 

V PPP+-/1  REDRNK  B ; []J(7  ; P.P  ; RB 

[1]  B SUBSTITUTES  POL., A,  FOR  THE  LAS'7'  VARIABLE  OF  B 

[2]  ->(1=  1 1 (pP  ) + QxRN+-pR3+(  -UIO+1  ) \pRES+B)  /OUT 

[3]  Ll  :-*■(  l<~l+pRES*-(  ( (PPpO  ) , “l  ) \ RES) ADD  A BY(RB  , ( 1 -~1  +p  RES) 

) + ( (P5«-“l+pPPP)  ,"1  )iRES)/Ll 

[4]  OUT:RES+(~llpRES)pRES 

V 
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RREDFAC [HIV 

V RES+-PLIB  REDFAC  RPOLY  ;FAC -,010  iTK  ;N  ;D  ;NFAC*  ;IF  ; 
Cl]  a ELIMINAmES  FACTORS  IN  PLIR  FROM  RPOLY (RAT . ”01 \ 

[2]  .w^pd  ,mK*-(i±pRPor,y)  )jt-ppoiY 

[3]  D+TKp (~1 ,TK)\RPOLY 
[41  NFACT++/ (l*p”LIB) 

[5]  RRL^-(ppPLIB)  -ni£«-l+IF«-0 

[6]  a SELECT  NEXT  FACTOR 

[ 7 ] NEWF  : -►(  NFACT<  IF*- IF+1 ) / mw  PF 

[ 8 ] FAC+-(  1 +PF4C  )p  F/4CVSIPIP  ( ( IF  -FF4  C"^  ) , FPL  p 0 ) + f ( IF - 1 ) 

PL  IB 

[9]  A IS  IT  A FACTOR  OF  NUMERATOR  AND  DENOMINATOR? 

[10]  £ : -►(  0 < + / , | (~l,(14-pf?))  +Q+-N  DIV  FAC)  /NEWF 

[111  -*-(0<  + /,  | ("1,(1  4pP))iP«-Z>  DIV  FAO/NEME 

[121  N+TKp(l ,TK+l\oQ)\Q 

[13]  D+TKp(l yTK+l\pP)\P 

[14]  +( (ppFAC)*+/pFAC)/L 

[15]  NOMORE:RES<-(PiN)  , [ 0 . 5 ] (P«-(ptf  ) Tpl) ) + D 

V 


P ; Q \RnTi 
) 


,RRLp  0 ) 4- 
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FILTER  PROBLEM  SOLUTION 

WN11  and  WD11  are  the  numerator  and  denominator  respectively  of  the  first 
diagonal  term  of  W,  expressed  as  a rational  polynomial  in  t and  0 . WT1  is  the 
same  term  expressed  as  a rational  polynomial  in  the  form  (2) , after 
elimination  of  common  factors  using  the  routine  REDFAC  and  the  library  of 
factors,  PLIB.  W1  is  the  same  term  after  elimination  of  6.  Similar  comments 
apply  to  the  other  diagonal  terms  of  W. 

RATPOL  takes  two  polynomials  as  input  and  returns  them  in  the  rational 

polynomial  form  (2),  after  cancelling  factors  which  are  powers  of  the 
variables 
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0 0 "4 

D+WT1+2 XPLIB  REDFAC  WN 11  RATPOL  WD  11 
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0 0 0 "6 
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0 0 0 0 
0 0 0 0 
2 "l  0 0 

2 0 0 0 
0 "l  0 0 

0 0 0 0 
0 0 0 0 

f>f/l<-(2_2pl  1 1 ~1)REDFAC(1  2 Ipl  1)REDRNK+WT1 
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Q+WT2+PLIB  RED FAC  WN 22  RAT POL  WD 22 

3 “2 
2 0 
1 2 
0 0 

2 “1 

4 "1 
2 “l 
0 "1 

Q+W2+(2  2 p 1 1 1 “l  )REDFAC(  1 2 ipl  1 ) RRMWK$Wm2 
1 0 "1 
14  1 

§WN3  3 


0 

0 0 

0 

0 

0 

0 

0 

0 

0 

0 

0 0 

4 

“12 

3 

8 

“12 

4 

0 

0 

0 0 
§WD  33 

“4 

24 

“6  0 

80 

“60 

2 4 

“4 

0 

0 0 

0 

"16 

0 

16 

“8 

0 

0 

0 0 

”4 

12 

“16 

16 

“1  2 

4 

5*PLIB  RED  FAC  WN  3 3 RATPOL  WD  3 3 

0 1 “l 

0 “1  4 

0 “1  “6 
0 1 4 

0 0 “l 

2 “1  0 
0 10 
2 “l  0 
0 10 
0 0 0 

□«.p/3«.(2  2 p 1 1 1 "1  )REDFAC(  1 2 Ipl  1 ) REDRNK<bWT3 
0 4 “4  0 

15  5 1 
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APPENDIX  III 

ORBIT  PROBLEM  SOLUTION 
C'Fll' ,[ 1 ] A) FORM  SEL  F 

5 4 3 4 3 5 

Fll=  10  2 3 V V - 878268  7 V + 93B60  7 VF  + 3 R7 6 1454  11  V 

3 3 3 2 2 7 

- 148  7 3 940  7 V W + 1189902  7 VF  - 318715236  U V 

2 5 2 3 2 2 3 

+ 234084492  11  V W - 53057340  U V W + 3479700  U VW 

9 7 5 2 

+ 654729075  UV  - 695674980  T)V  11  + 266431410  UV  V 

3 3 4 

- 42723300  UV  W + 2320275  777 


5 4 2 4 3 4 

711=  - V + 53640  U V -1009  V W - 6379326  77 

32  32  2 fi  24 

+ 1657140  7 V F - 40446  U F + 97257888  U V - 57948120  7 V F 

2 2 2 2 3 9 p 

+ 9297840  7 V F - 255000  7 .¥  - 310194825  UV  + 290768940  UV  V 
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Figure  1.  Formation  of  the  product  of  polynomials 

of  a single  variable  from  the  outer  product 


Figure  2.  Formation  of  polynomial  product  by  rotation 
of  outer  product.  The  padded  outer  product 
is  skewed  along  its  last  axis,  and  then 
reduced  by  summation  along  its  first  axis. 
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